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Abstract 

On the base of the diffusion Monte-Carlo method we develop the method 
allowing to simulate the quantum systems with complex wave function. The 
method is exact and there are no approximations on the simulations of the 
module and the phase of the system's wave function. In our method averaged 
value of any quantity have no direct contribution from the phase of distri- 
bution function but only from the phase of the Green function of diffusion 
equation. 

We test the method by the simulations of the ground state of fermions in 
two-dimensional parabolic well. Anyons are used for the representation of the 
two-dimensional (2D) fermions. We vary the number of fermions from two to 
ten and hnd a good agreement of the numerical results with analytical ones 
for the numbers of the particles N > 4. 

I. INTRODUCTION 



There are several objects of the quantum mechanics where the investigated system's wave 
function is essentially complex (having imaginary part). The examples of these objects are 
a two-dimensional system of electrons in the external uniform magnetic field when vector 
potential has a central-symmetric form and system of anyons. The complexity of wave 
function does not give a possibility for the simulations of such kind of systems by using 
well-known Green Function Monte-Carlo method (see review This method essentially 
demand the reality of the system's wave function which considered as probability weight in 
during stochastic process. 

There had been undertaken the several attempts [Q-|J for the construction a method for 
the simulation of quantum systems with complex wave function. Authors of the works |]2]|| 
suggested to take as a probability weight the module of the complex distribution function. 
All quantities are calculated in |2|U by averaging over this complex distribution function. 



1 



The phase of this function 7 is taken accordingly 7 = where 7« is the phase at the 

i-time step of the stochastic process. 

Another quantum Monte Carlo method using algorithm without branching for the sim- 
ulating of complex problems was developed in [|J. The main trouble of this approach was 
an increasing of the statistical error as a function of the whole time of the simulation. 

In the previous paper || it was suggested new Complex Diffusion Monte-Carlo (CDMC) 
method allowing to simulate the quantum systems with complex wave function. Our method 
is rather close to the fixed - phase Diffusion Monte-Carlo method developed for the simulation 
of the two-dimensional electrons in the magnetic field 0. In this work there had been an 
exact simulation only of the module of system's wave function. The phase of the wave 
function was considered as fixed and equal to phase of the Laughlin's wave function [0. In 
contrast, our method include also the simulation of the phase factor of the wave function. 

The basic difficulty of the numerical simulation of the fermions is essentially the same as 
for the systems with the complex wave function. Their wave function can change the sign 
and therefore can not be used as the probability weight in the simulation process. 

For the simulation of the continuum (not on the lattice) fermionic systems it were de- 
veloped widely used fixed node Monte-Carlo method ||( see also reviews |]||) and recently 
proposed constrained path Monte-Carlo method [ID]. In these both methods it was assumed 



the restriction on the random walks connected with the uncertainties in the space localiza- 
tion of the wave function node surfaces. The comparison of these methods was done in 
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In the 2D space we have unique tool - so-called anyons, particles (bosons or fermions) 
with additional gauge interaction which provide needed statistical properties [12||13|] . So, by 
tuning the coupling constant we may arrive to the fermions starting from the bosons. It looks 
very attractive to apply the anyons for the simulation of the 2D fermions, because statistical 
vector potential that gives anyonic property is relatively smooth and continuous []T3| , |T^| . The 
remaining main problem is that the wave function of anyonic system is essentially complex, 
i.e. it contains an imaginary part. 

The CDMC method was successfully tested by the simulation of ground state of one 
electron in magnetic field || . In the present work we briefly outline the idea of CDMC (the 
details of the algorithm one can find in ||) and provide a further test it by simulating the 
ground state of fermions in 2D parabolic well and comparing with well known analytical 
answer for the ground state energy as a function of the number of fermions. The simulation 
is done for the fermionic systems with the number of particles from two to ten. It was 
found a good agreement of the numerical results with analytical ones for the numbers of the 
particles great than four. This observation of the fermionic simulation allows us to hope 
that the CDMC is a good tool for the simulation of 2D fermion systems. 

II. MODEL AND CDMC METHOD FOR THE SIMULATION 

Hamiltonian of anyons in the 2D parabolic well has a form: 

1 N - Mw 2 N 

H=^7 £(* + mf + ^ E rl (1) 

i=l i=l 
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Here M is the mass of particle, p = — zftV, where V = w Q stands for characteristic 

ox ay 

frequency of free particles in 2D parabolic well and is radius vector of z'-th particle. The 
number of particles is N. 

Vector potential for anyons A(fi) [13|,[14| in (1) is 

N -> -> 

m) = nuj: z -^f- (2) 

j>i 

Here z is unit vector perpendicular to ID plane and v characterizes the form of fractional 
statistics ( spin of the anyon ). In this work we want to employ a bosonic representation of 
anyons, therefore v = corresponds to noninteracting bosons and v = 1 is the case of our 
interest - free fermions. 

In CDMC a Monte Carlo simulation is carried out with the complex distribution function 

f(R,t) = ** T (RMR,t), (3) 

where R stands for the coordinates of all particles, ^t{R) is a trial wave function and the 
wave function ty(R,t) satisfies a Schrodinger equation with imaginary time (expressed in h 
units) 

-^^ = (ti-E T MR,t). (4) 

With a suitable choice of a trial energy Et it is possible to arrange the convergence 
*&(R, t) — > ^o(R) - an exact stationary wave function of ground state of Hamiltonian ([[]) at 

For a bosonic representation of anyons we take (conjugated) trial wave function in the 
form: 

N 

n(R) = um?i)- (5) 
i=i 

Good choice of single particle trial wave function ^f T (fi) ( see below ) as energy Et must 
provide a convergence of simulation process to exact ground state. 
The distribution function / = f(R,t) satisfies a diffusion equation: 



N N 

- = -D 

dt 



° j -D £ A,/ + D £ V^/ReF^rO) + 

i=l i=l 



N - - h - - 

+z£[V i p/IrnF Q (r i )) - ^A(r i )V i f] + (E L (R) - E T )f. (6) 

i=l 

h 2 

Here D = tttJ; = . When a time step integration of diffusion equation (^), r, goes to 
zero, then Green function for this equation has a form: 
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G(R,R';t) 



explDr^AKfi,^ ')} 

<=i x 

(AttDt) n 

N 

^(rl-f/- J DrReF Q (r i ')) 2 



x exp 



x exp 



i=i 



ADt 



exp 



-t{E l {R)-Et) 



x 



N 



]T ^(n, r, ')(r< ~n'- DrReF^ ')) 



8=1 



In the expressions (|6|) and (^) 



F g (f i ) = 2**- i (r i )V i ^(n), 



(7) 



(8) 



£ L (#) = fj" (iJ)ff $r(iJ), 



(9) 



2 JV 



E(p*-^)) a + - 5 - £ E*f- 



ZiW i=l z i=i 



In general case ^(fi) is a complex wave function and 

F Q (fi) = ReF Q (r i )+iImF Q (r i ). 
In flj]) (see also 0) we have introduced a new quantity 

A Q (n^i = 5lmF Q (? ; ') - gA(?i). 



(10) 



(11) 



(12) 



In the expressions (|7|) and ( |I2| ) vectors i? and r*j correspond to time point t + r and vectors 
R f ' to time point £. 

One can see that the Green function G(R, R '; r) ([/]) as and distribution function / in 
(0) is a complex. These both quantities are satisfying to the usual integral relation 



f(R, t + r)= dR 'G(R, R r)f(R t). 



(13) 



From relation (|13 ) it is following that the module and phase of the distribution function at 
consequent time point are determining by the module and the phase of the Green function 
and of ones of the distribution function at previous time point of integration of diffusion 
equation (H). 

— * — * 

As the energy E L (R) is a complex, real and imaginary part of E L (R) contribute separately 
to ones of the Green function. Therefore the expression of the last two exponents in (|7]) has 
a form: 



exp 



x exp 



N 



-r(ReE L (R) -E T ) 



x 



]T A Q {r h n ')(n - r, ' - J DrReF Q (r i ')) - irImE L (R) 



i=l 



(14) 
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The mean value of the some quantity B(R(t)) on the complex distribution function 
f(R,t) at the time t of the simulating process is calculated by the formula (see P||) 



M 



^exp[z 7 ( J R 4 (t)] J B( J R J (t)) 

< B{t) >= ^-^ . (15) 

f>xp[z 7 (%)] 
i=i 

Here Ri(t) is the coordinates of all particles from the z-th configuration of the M - config- 
urations at time t. The coordinates Ri(t) of particles in (|15|) are weighted with probability 



\f(Ri,t)\ and the quantity j(Ri(t)) is a phase of the distribution function f(Ri,t). 

From integral expression ( fL3|) it is seen that phase of the distribution function •y(R i (t+T)) 
connected with the phase of the Green function ja(Ri(t + t), Rj(t)) and the phase ^(Rj^t)) 
of the distribution function by the relation: 

l(Ri(t + r)) = lG (R t (t + T),Rj(t)) + 7(£i(t))- (16) 

Here Rj (t) is j- th configuration at time t. So, we can write 

M 

y^ e *7c?(^(t+r),-^(t))+i7(4-(*)) B(Ri(t + t)) 

<B{t + r)>=^ 



M 



i=i 

2 e i ^^ t+T ^^B{R i {t + r)) 



i=i 

M 



(17) 



i=i 



M 

y^ e i7G(^i(t+r),^jW) 
t=l 

The expression (^) is hold, because all quantities in numerator and denominator are 
weighted with probability \f(Ri(t + r),t + t)\, i.e. at time moment t + r with new configu- 
rations -Rj(t + r). Hence, a mean quantity < B{t + r) > is determined only by phase 7g of 
the Green function. It is a key conclusion, providing the calculation of mean quantities in 
our CDMC. 

Next, by taking into account the expression for the simulation of the random displacement 
for the particle r, \ = r, \ ' + DrReFQ^ ') + 2(Dr) 1 / 2 % (see and ||), where x is a gaussian 



random number, the expressions for Fq(fi) @, for Ei(R) (|9]) and for Aq^, ') ([121) and 
also the expression ([14]), one can show that at r — >• the phase of the Green function jo 
tends to zero as r 1 / 2 . This property of 7g together with the formula ([T7|) allows CDMC 
successfully simulate all mean quantities. 

The simulation process starts with mean field approximation of Fetter, Hanna and Laugh- 

lin [^] . In this approximation the mean field A generated by the average density p of the 
particles: 

A(f) = pirhv(zx r) = ^B x r. (18) 
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Here B = l-nphvz is fictitious uniform mean magnetic field that defines corresponding 



magnetic length an = yh/B and cyclotron frequency oo c = B/M. Particles move in this 

mean magnetic field. Energy spectrum of Hamiltonian ([]]) with A and without 2D parabolic 
well is Landau levels. We have to simulate anyonic ground state. Therefore we take a wave 
function ^(rj) in the following form: 

•ko-M-^M-^). (19) 

Here C is a normalization constant, a is numerical parameter and R Q = (h/Mwo) 1 ^ 2 . In 
this expression first exponential function with a = 1 corresponds to the ground state wave 
function of particle in 2D parabolic well and second one to the ground state of Landau 
levels. 

Average density of particles p is given by p = l/rcr 2 where r Q is the mean distance 
between particles in system. We assume that r Q equals to R a . The reason of this is that the 
interaction between two anyons at any non-zero v has repulsive character and the area S n 
that occupied by free particle at n-th 2D oscillator's level is S n = mrR 2 . Last relation for 
S n follows by using virial theorem for the quantum energetic states for the free particle in 
2D parabolic well. 

By substituting expression B into expression a# and taking account that p = 1/rcr 2 , we 



find the following relation between a# and r : a# = ■Jl/2ur t} . By similar way one can show 
that hu c /2 = vhu . 

It is natural to use hu c /2 as energy unit, an as length unit and 2/ftui c as time unit r. In 
these units the last term of the Hamiltonian flU) has a coefficient l/4z/ 2 and the argument 
of the first exponential function in (|l!|) has a coefficient a/Av. By using above relations 
between %uj c /2 and Tiujq and between an and ro it is easily to express the energies in terms 
of 2D oscillator quanta hu>o and the lengths in terms of tq. 



III. DISCUSSION OF THE RESULTS 

The simulation of the ground state of the systems of several fermions in 2D parabolic 
well was performed. By setting the coupling v — 1 in (0) we get fermions from anyons. This 
value of v is fixed below. As in || , the total simulation time t is divided into several (in our 
case we took 40) time blocks At. In every time block At initial number of configurations 
N c is chosen equal to 1000 and number of time steps r equal to 30. 

We simulate the energy Esim, mean radius r and mean square of radius of r 2 of the 
fermionic systems. The calculation of these quantities are performed by the formula (|l~7|). 
Here we take into account that the energy E L (R) has as a real part ReE L (R), and also an 
imaginary part IulEl(R). The number of fermions is varied from two to ten. 

The t dependence of the energy E S im has a general structure A + Ait 1 / 2 where A 
and A\ are the numerical constants. 

We have obtained that at a very small r the mean deviation of Esim larger than the 
fluctuations of Esim- At large r further increasing of it leads to the increasing of the 
fluctuation of E$im an d thus to the growing of the population number N p (the population 
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number N p is connected with first exponential function in (|T4] ) and was also introduced in 
H). The very large values of r give a dramatical increasing of N p that leads to the instability 
of E SIM and whole stochastic process. So, we have to choose optimal time step r . 

CDMC has an artifact that the averaged real quantities become complex, but real 
parts of them are essentially bigger than imaginary parts because the imaginary parts are 
controlled by the phase of the Green function jo depending as Tq 5 . Then the imaginary 
part of Esim for almost all number of fermions iV is less than 0.15 percent of its real part. 
For f and r 2 they are less than 0.07 percent of their real parts. 

Table 1 present the exact analytical result for the ground state energy of the fermionic 
systems in the 2D parabolic well E EX act, E S im (in ftw units) and the simulated values of r 
and r 2 ( in r Q and r 2 units respectively). The values of E EX act was calculated by occupying 
the available lowest states. The energy of K-th quantum number state of the 2D oscillator 
is Ek = hco (K + 1) and degeneracy number of this state is equal (K + 1) fill . 

Esim is depended on value of numerical parameter a in the trial wave function ty^i^i) 
(|19|) . We take a providing minimal value of the energy Esim- The numerical data of a 
that give a minimal energy Esim are outlined also in the Table 1. We have compared the 
CDMC dependence of the parameter a versus number of particles iV with the variational 
calculation of a f!7 |. If at the changing of N from two to ten CDMC a changes from -0.593 
to -0.084, so variational a changes from -0.293 to -0.683 (see flT]]). For the big number 
fermions (as one can see from table) CDMC a is almost the same for all systems. Thus a 
mean field approximation of Fetter, Hanna and Laghlin |15| with numerical parameter in 
a wave function can be considered as a good start approximation for the simulation of a 
big number anyons. Table 1 also shows the numerical values of the population number N p 
and optimal time steps r . Good agreement of the simulation results for the Esim with the 
exact ground state energies E EX act for the numbers of the particles bigger than four allows 
us to hope that CDMC gives a possibility for the exact simulation of the big number 2D 
fermions. 
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TABLES 



TABLE I. The energies of the fermion systems in 2D parabolic well. 
Here iV - numbers of particles, Eexact - analytically calculated ground state energies (in hu> 
units), Esim - results of the simulation ground state energies (in %uj units), r - the numerical 
values for the simulated mean radius (in r a units), r 2 - the numerical values for the simulated 
mean square radius (in r 2 units), a - the numerical parameters in the trial wave function (19) 
that gives a minimum Esim, N p - the population numbers and tq - the optimal time steps. All 
simulated quantities and their deviations from mean values were averaged over 30 last time blocks 
At when Esim and N p have had relatively stable values. 



N 


Eexact 


Esim 


r 




a 


N p 


TO 


2 


3 


1.842 
±1.297 • 10~ 3 


0.935 
±6.994 • 1(T 4 


0.8852 
±1.319 • 10" 3 


-0.593 


980 ±5 


0.08 


3 


5 


4.399 
±4.471 • 10~ 3 


1.092 
±1.564 • 10" 3 


1.573 
±3.814 • 10" 3 


-0.423 


925 ±9 


0.06 


4 


8 


7.399 
±4.092 • 10~ 2 


1.236 
±1.844 • 10" 2 


1.919 

±4.562 • 10" 2 


- 0.3 


908 ± 13 


0.04 


5 


11 


11.195 
±1.798 • 10~ 2 


1.452 
±3.738 • 10~ 3 


2.614 
±1.324 • 10" 2 


-0.453 


949 ± 10 


0.04 


6 


14 


14.492 
±3.136 • 10~ 2 


1.619 
±4.261 • 10~ 3 


3.275 
±1.564 • 10" 2 


-0.092 


992 ± 24 


0.02 


7 


18 


18.446 
±9.126 • 10~ 2 


1.752 
±8.277 ■ 10" 3 


3.855 
±3.665 • 10- 2 


-0.022 


1042 ± 53 


0.02 


8 


22 


22.383 
±8.569 • 10~ 2 


2.020 
±5.821 • 10" 3 


5.112 
±2.790 • 10" 2 


-0.046 


1110 ±57 


0.015 


9 


26 


26.556 
±9.595 • 10~ 2 


2.294 
±1.061 • 10" 2 


6.672 
±6.435 • 1(T 2 


-0.067 


1090 ± 33 


0.01 


10 


30 


30.905 
±2.010 • 1Q- 1 


2.601 
±2.839 • 1(T 2 


8.599 
±2.095 • KT 1 


-0.084 


1032 ± 23 


0.005 
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